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ABSTRACT 

The measure of the third-order structure function, T, is employed in the solar wind to compute the cascade 
rate of turbulence. In the absence of a mean field Bo = 0, T is expected to be isotropic (radial) and independent 
of the direction of increments, so its measure yields directly the cascade rate. For turbulence with mean field, 
as in the solar wind, Y is expected to become more two dimensional (2D), that is, to have larger perpendicular 
components, losing the above simple symmetry. To get the cascade rate one should compute the flux of Y, 
which is not feasible with single-spacecraft data, thus measurements rely on assumptions about the unknown 
symmetry. We use direct numerical simulations (DNS) of magneto-hydrodynamic (MHD) turbulence to char¬ 
acterize the anisotropy of Y. We find that for strong guide field B 0 = 5 the degree of two-dimensionalization 
depends on the relative importance of shear and pseudo polarizations (the two components of an Alfven mode 
in incompressible MHD). The anisotropy also shows up in the inertial range. The more Y is 2D, the more 
the inertial range extent differs along parallel and perpendicular directions. We finally test the two methods 
employed in observations and find that the so-obtained cascade rate may depend on the angle between Bo and 
the direction of increments. Both methods yield a vanishing cascade rate along the parallel direction, contrary 
to observations, suggesting a weaker anisotropy of solar wind turbulence compared to our DNS. This could be 
due to a weaker mean field and/or to solar wind expansion. 

Subject headings: The Sun, Solar wind. Magneto-hydrodynamics (MHD), Turbulence. 


1. INTRODUCTION. 

Magneto-hydrodynamic (MHD) turbulence in the presence 
of a mean-field Bo has a tendency to become two-dimensional 
(2D). This tendency was early recognized by inspection of 
the Fourier energy spectra in direct numerical simulations 
(DNS). The energy distribution is indeed anisotropic, resid¬ 
ing in wavevector mostly perpendicular to the mean mag¬ 
netic field (Montgomery & Turner 1981; Shebalin et al. 1983; 
Grappin 1986). Ideally one would like to quantify the two- 
dimensionalization as a scaling relation between parallel and 
perpendicular wavenumbers having the same energy density, 
fell o' kj. If p - 1 the anisotropy is scale independent, and the 
aspect ratio fey/fe ± of the isocontour of the Fourier spectrum 
does not change with scales. If instead p < 1 the aspect ra¬ 
tio increases with wavenumber, that is, the spectrum becomes 
more and more 2D at smaller and smaller scales. In DNS, 
the parallel spectral extent is generally very short, due the 
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limited achievable Reynolds numbers, and the parallel spec¬ 
trum rarely shows a power-law, rendering the distinction be¬ 
tween scale-dependent and scale-independent anisotropy dif¬ 
ficult. In contrast, the two-point correlation in real space ex¬ 
pressed by II-order structure functions, S , shows in general 
nicer power-law scaling in the parallel direction, allowing one 
to quantify the scale-by-scale anisotropy. In analogy with 
Fourier spectra, the scaling relation involves parallel and per¬ 
pendicular increments that have the same energy l\\ oc t'j, also 
known as eddy anisotropy. 

The II-order structure function anisotropy has been widely 
studied in DNS of incompressible MHD turbulence. Us¬ 
ing a local mean-field to identify parallel and perpendicular 
increments, one finds a scale-dependent anisotropy (Cho & 
Vishniac 2000): the anisotropy grows at smaller and smaller 
scales, suggesting a complete 2D state at small enough 
scales. Furthermore, the anisotropy is controlled by b rms /Bo 
where b rms indicates the root-mean-square amplitude of turbu¬ 
lent fluctuations, the stronger Bo the stronger the anisotropy 
(Muller et al. 2003). The anisotropy is also stronger in re¬ 
gions of stronger magnetic field (Milano et al. 2001). Sim- 
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ilarly, the anisotropy in the Fourier spectrum increases for 
the stronger Bo (Oughton et al. 1998). However, employing 
structure function to measure the anisotropy with respect to 
the global mean-field returns a scale-independent anisotropy 
(Chen et al. 2011), implying that the two-dimensionalization 
does not increase at smaller scales but reaches an asymptotic 
value. A similar dichotomy exists in solar wind measure¬ 
ments, in which one can compute the two-point correlation 
in time from time series of data collected in-situ by space¬ 
craft and then adopt the Taylor hypothesis to obtain spatial 
increments. The increments are thus taken along the radial di¬ 
rection, but the anisotropy with respect to the magnetic field 
is recovered thanks to its variable direction with respect to 
the radial. As in DNS, the structure function S is found to 
be more energetic along perpendicular increments than along 
parallel increments. Again, a local mean-field analysis yields 
a scale-dependent anisotropy (Horbury et al. 2008; Wicks 
et al. 2010; Chen et al. 2012), while a global mean-field anal¬ 
ysis indicates a scale-independent anisotropy (Tessein et al. 
2009). Several authors (Cho et al. 2002; Chen et al. 2011; 
Beresnyak 2012) showed that for strong turbulence the scale- 
dependent anisotropy is smoothed out in a global mean-field 
analysis, even in the presence of a strong mean-field. On the 
other hand, Matthaeus et al. (2012) noted that since the lo¬ 
cal mean-field is a random variable, the local II-order struc¬ 
ture functions involve higher-order statistics and cannot be the 
real space equivalent of the power spectra. However for small 
enough b rms /Bo the global and local measures are expected to 
coincide. 

In this work we investigate the process of two- 
dimensionalization of MHD turbulence, focusing on the 
anisotropy measured in the global frame. In this frame, one 
can obtain a dynamical equation (labelled KHYPP equation) 
that relates II-order and Ill-order structure functions (Politano 
& Pouquet 1998), extending to incompressible MHD the Von 
Karman-Howart-Monin equation for incompressible hydro- 
dynamic turbulence. According to the KHYPP equation, for 
stationary and homogenous turbulence in the inertial range, 
the divergence of the Ill-order structure function Y is propor¬ 
tional to the cascade rate of turbulence 4e = -V • Y. The 
divergence is negative, implying that the cascade is achieved 
by removing positive correlations, and thus increasing the 
amplitude of S (flattening its power-law index). Thus, by 
characterizing the anisotropy of Y one can get insight into 
the process of two-dimensionalization of MHD turbulence. 
Previous studies on the Ill-order structure functions in DNS 
of the MHD equations were limited to 2D (Politano et al. 
1998; Sorriso-Valvo et al. 2002) thus leaving out the issue 
of anisotropy. In a recent work, Lamriben et al. (2011) re¬ 
ported for the first time the vector Y measured in an experi¬ 
ment of rotating hydrodynamic turbulence. As rotation was 
increased, the anisotropy of the II-order structure functions 
also increased. They found that the two-dimensionalization 
can be associated with the tilting of the vector Y toward the 
plane orthogonal to the rotation axis and that the tilting begins 
at small scales and then propagates to larger and larger scales. 

In the present work, we carry out a similar analysis on data 
from three-dimensional (3D) DNS of incompressible MHD 
turbulence by computing for the first time the 3D Ill-order 
structure functions, in the presence or absence of a mean-field. 
We find that the degree of two-dimensionalization as mea¬ 
sured by S is associated to the relative excitation of pseudo 
and shear Alfven polarizations for stationary turbulence with 
mean field Bq. We also analyze the full KHYPP equation in 


developed decaying isotropic turbulence, offering a descrip¬ 
tion of the cascade in real space based on structure functions 
(SF), analogous to the usual Kolmogorov cascade in Fourier 
space. Note that while the latter is based on the assumption 
of locality, the KHYPP equation is free from this assumption 
(requiring only homogeneity), thus representing a more gen¬ 
eral description of the cascade process in MHD turbulence. 

The Ill-order structure function Y has also been computed 
in solar wind data to obtain the cascade rate of solar wind 
turbulence (Sorriso-Valvo et al. 2007; Marino et al. 2008, 
2012; Macbride et al. 2008; Smith et al. 2009; Stawarz et al. 
2009, 2010). These rates are consistent with the heating rate 
estimated from proton temperature gradients (Vasquez et al. 
2007; Cranmer et al. 2009), suggesting that turbulence may 
supply the heating required to sustain the non-adiabatic ex¬ 
pansion of the solar wind. However, applying the KHYPP 
equation to the solar wind is a bit problematic since solar 
wind turbulence is neither stationary nor homogenous (e.g. 
Hellinger et al. 2013; Gogoberidze et al. 2013; Dong et al. 
2014). To obtain the cascade rate one should compute the di¬ 
vergence of y, which is quite difficult with single-spacecraft 
data, since increments are taken only along one direction at 
time (see however Osman et al. 2011 for an integral form that 
exploits the four CLUSTER spacecraft with the minimal as¬ 
sumption of axisymmetry). The cascade rate can thus be re¬ 
trieved only assuming a form for the unknown anisotropy of 
y. Although some theoretical predictions exist (e.g. Podesta 
et al. 2007; Galtier 2009), two methods are commonly em¬ 
ployed in observations that assume respectively isotropy or 
an anisotropic model based on the geometrical slab-plus-2D 
turbulence that was introduced by Matthaeus et al. (1990) 
to describe the two-point correlation function of solar wind 
turbulence. We will exploit the data from our DNS to test 
the two methods employed in solar wind data against known 
anisotropic Ill-order structure functions and estimate the pos¬ 
sible systematic errors. 

The plan of the paper is as follows. In section 2 we give a 
brief introduction to the KHYPP equation, while in section 3 
we describe the method employed to compute structure func¬ 
tions (SF) of II-order and Ill-order. The results are presented 
in section 4, where we first describe the anisotropy of II-order 
structure functions of the simulations considered. The rest of 
the section is dedicated to the Ill-order SF. We consider first 
a simulation of decaying turbulence without mean magnetic 
field, allowing us to test the soundness of our analysis method 
and to verify that the time-dependent KHYPP equation holds 
in developed decaying turbulence. Then, we analyze two sim¬ 
ulations of turbulence with mean-field that have a different 
strength of anisotropy. Finally, in section 5, we test on runs 
with Bo + 0 the two methods employed to measure the cas¬ 
cade rate in the solar wind. We conclude with a discussion on 
the results and on the application to the solar wind turbulence. 

2. STRUCTURE FUNCTIONS AND THE KHYPP 
EQUATION. 

The Von Karman-Howart-Yaglom, Politano-Pouquet equa¬ 
tion (KHYPP) for non-stationary, anisotropic, and incom¬ 
pressible MHD (Politano & Pouquet 1998; Podesta 2008; 
Carbone et al. 2009) can be obtained from the original MHD 
equations written in term of the Elsasser variables z ± = 
u + b/ oj4np, by subtracting the MHD equations evaluated at 
different positions x and x + C and by averaging in the volume. 
Under the assumptions of incompressible, homogeneous tur- 
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bulence one obtains 


5,<|Az ± | 2 > + V f -(Az T |Az ± | 2 ) = 
-II - A + 2vV»(|Az ± | 2 ) - 4e ± , 


( 1 ) 


where we have defined the two-point correlation, A z ± (x, t ) = 
z ± (x + {) - Z ± (x), and (...) stands for the volume average. 
These equations describe the evolution of the II-order struc¬ 
ture function for each Elsasser variable, S ± = <|Az* | 2 ). The 
divergence term in the left hand side is the Ill-order struc¬ 
ture function, Y ± = ^Az + |Az ± | 2 ^, involving products of Az + , 
and Az _ , which we name Yaglorn flux in the following. On 
the right hand side (rhs), n and A represent pressure terms 
and sweeping terms (responsible for the Alfven effect) respec¬ 
tively, both vanishing for globally homogeneous turbulence 
(Carbone et al. 2009). The remaining terms represent dissi¬ 
pation. The first one involves the Laplacian with respect to 
the increments t, it vanishes for vanishing viscosity (for sim¬ 
plicity we assumed equal viscosity and resistivity v = rj). The 
second one, e ± = -d,E ± = v(T.j(djzf)(djzf)), is the dissi¬ 
pation rate of the Elsasser energies (E ± = (|z ± | 2 /2)). In the 
former, the derivatives of the primitive fields zr do not com¬ 
mutate with the averaging operation, and the dissipation rate 
remains finite for vanishing viscosity. 

Summing the contributions of both Elsasser fields, one fi¬ 
nally gets an expression for the total energy and cascade: 


d,S + V € ■ Y = —4e + 2vV 2 S, 


( 2 ) 


where 5 = 1/2(5 + + S~), Y = 1/2 (Y + + Y~), and e = 
1 /2(e + + e~). Note that because of homogeneity in the above 
Equations (l)-(2) all the variables depend only on the vector 
separation, t. 

This equation is valid for decaying turbulence and describes 
the classical scenario of a turbulent flow in which the dissipa¬ 
tion of energy is achieved through a cascade of energy toward 
smaller scales, where fluctuations are finally damped by vis¬ 
cosity. In forced turbulence one should add on the right hand 
side the forcing terms {J r ) that inject energy (usually) at large 
scales. 

For stationary turbulence ( d,S — 0) forced at large scales 
CF ^ 0 only at large scales) the injection, cascade, and dis¬ 
sipation all occur at the same rate. At high Reynolds number 
one expects their respective ranges to be well separated, in 
analogy to the Kolmogorov cascade in Fourier space, so one 
has that: i) at large scales the second and last terms in Equa¬ 
tion (2) are negligible and T = 4e, the forcing balances the 
dissipation rate, ii) at small scales the second term is negligi¬ 
ble, 2vV 2 5 = 4e, and the damping rate is equal to the dissipa¬ 
tion rate, and finally iii) at intermediate scales the last term is 
negligible, yielding 


Vf Y = —4e, 


(3) 


that is, the cascade rate, which is equal to the dissipation rate, 
is given by the divergence of the Yaglom flux. Note that Equa¬ 
tion (3) can be used as a definition of the inertial range as 
being the ensemble of scales for which the equation is ap¬ 
proximately satisfied. The definition should hold for quasi¬ 
stationary forced turbulence and for developed decaying tur¬ 
bulence. As we will see, for developed decaying turbulence, 
the time-dependent term is non-negligible only at large scales 
where d,S = -4e (the dissipation is balanced by the decay of 
the II-order SF at large scales). 


We can now give a more physical interpretation of the cas¬ 
cade process by rewriting Equation (2) in terms of the autocor¬ 
relation function, C = C'+C .with CAT) = (z ± (x+()-z ± (x)). 
The autocorrelation functions are related to structure func¬ 
tions by 

S ± (() = 2E ± -2C ± ((), (4) 

and using <3, E ± = —e ± , the KHYPP equation becomes: 

d,C - V e ■ Y = -2vV 2 C, (5) 

showing that Y is a flux of negative correlations. A permanent 
flux of negative correlations toward small scales is equiva¬ 
lent to constantly building new small-scales gradients. For 
instance, the formation of 2D quasi-perpendicular turbulence 
will be revealed by a Yaglom flux bringing negative correla¬ 
tion at small perpendicular scales, hence the vector Y must be 
quasi-uniform, parallel to the i ± axis, and pointing toward the 
parallel €\\ axis. 

Coming back to Equation (3), for turbulence in stationary 
state or decaying turbulence, the cascade is a constant at all 
inertial-range scales (-V • Y = const). Thus, the inertial 
range anisotropy cannot appear as different cascade rates in 
the parallel and perpendicular directions. The anisotropy in¬ 
stead will show up in the shape of the domain of ( for which 
V • Y = const. 

To illustrate such an anisotropy, one can assume some par¬ 
ticular symmetry of the flux Y to characterize the cascade, 
with the additional advantage of obtaining a direct relation to 
the cascade rate, so as to avoid computing the divergence of 
Y. The simplest assumption is that of isotropic turbulence 
for which Y depends only on the scalar increment t. Rewrit¬ 
ing the divergence in spherical coordinates, and assuming sta¬ 
tionary conditions and vanishing viscosity. Equation (3) in the 
inertial range becomes the isotropic KHYPP equation, yield¬ 
ing 


j«, 3 Y_m 

L 4 l ’ 


(6) 


in which the Yaglom flux Y( = Y ■ l/\t\ is projected along the 
increment. This form is often used in solar wind studies, since 
although one does not have access to the full divergence in in- 
situ data, the cascade rate can be obtained directly from the 
projected Yaglom flux. The inertial range occupies a volume 
that is a sector of a sphere, it can be defined as isotropic since 
it has the same extent and location on parallel and perpendic¬ 
ular increments. 

A strongly anisotropic case case is that of 2D turbulence, 
obtained when the Yaglom flux in the inertial range depends 
only on the in-plane increments, t ± : 


V Y = V x Y ± = —4e, 


(7) 


where Y ± are the in-plane components of V and V ± denotes 
derivatives with respect to the in-plane increments. The Ya¬ 
glom flux can have out-of-plane components, but the cascade 
rate is determined only by the in-plane components. Note that 
a Yaglom flux having only in-plane components in the inertial 
range is a sufficient condition to have a 2D cascade. Assum¬ 
ing isotropy of the in-plane increments, that is, a dependence 
only on the scalar separation € ± , one obtains again a direct 
relation between Ill-order SF and the cascade rate: 


e 2fl = - 


YeAfJ 

2( ± ’ 


( 8 ) 
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TABLE 1 

Runs and parameters for simulations. 

forcing 
decaying 

k jf = 1/5 
< 2 
|V) < 2 

frozen 

Note. — b rms = is the root mean square magnetic field fluctuation. 
Bo is the mean magnetic field along the x axis. R x = L x /Ly is the aspect ratio 

of the box of size L y = L z = 2 n. The parameter^ = t^// NL = l/ ± b rms /k^Bo 
controls the strength of turbulence at forcing scales. N x , N y , and N z are the 
number of grid points, v is the viscosity coefficient (equal to the resistivity 
77 ). Re = [2n /(k^Ldiss)] 4 ^ is the effective Reynolds number, where the dis¬ 
sipation scale is defined as Ldi SS = (v 3 / 6 )^ 4 . For run A, Re decreases with 
time in the indicated interval. Forced wavenumbers are normalized by L y , 
with k\\ = k x , k± = ( k y + k z ) 1 ^ 2 . In run B the sound speed is Cs ~ 12 and the 
conductivity coefficient is k = v. 


Run brms Bo R x X N x ■ N y ■ N z 10 W Re 
A 05 0 T - 1024 3 09 2.6 

-» 2.2 

B 09 5 5 1.8 512 3 1.5 1.1 

C 1 5 1 0.2 256 ■ 1024 2 1 2.8 


with Y( ± — Y ■ (J\( ± \ indicating the projection of Y ± on 
the radial direction in cylindrical coordinates. Note that the 
inertial-range domain is not confined to the 2D plane even in 
this anisotropic turbulence. As we will see, the anisotropy of 
the inertial-range domain shows up in its different extent and 
location along parallel and perpendicular increments. 

For completeness, we consider finally the case of ID turbu¬ 
lence, when the Ill-order SF depends only on one coordinate, 
say, £\\. One obtains the cascade rate as 


id _ 

" “ 4f„ 


(9) 


Note that the geometrical model of slab turbulence does not 
correspond to the ID turbulence, since in the former, fluc¬ 
tuations are assumed to be perpendicular but to depend only 
on parallel wavevectors (and hence parallel separations). The 
slab geometrical configuration would indeed have a vanishing 
divergence. 


3. SIMULATIONS AND NUMERICAL METHOD. 

We consider three high-resolution simulations of MHD tur¬ 
bulence whose parameters are listed in table 1. Run A is a de¬ 
veloped decaying simulation of incompressible MHD turbu¬ 
lence without mean-field, representing isotropic turbulence. 
Run B is a simulation of weakly compressible MHD turbu¬ 
lence 2 , with mean-field Bq — 5 and anisotropic forcing. The 
forcing is applied only to components perpendicular to the 
mean-field and to wavevectors mainly perpendicular to the 
mean-field. Thus we are dealing with strong anisotropic tur¬ 
bulence of fluctuations with shear-Alfven polarizations, a con¬ 
figuration akin to Reduced MHD 3 . Finally run C is a simula¬ 
tion of incompressible MHD, again with mean-field Bq — 5, 
but with a forcing which is isotropic in both components and 
wavevectors. In the latter the forcing is actually a freezing 
of the modes 1 < k < 2 that maintains an equal amount of 
pseudo-Alfven and shear-Alfven polarizations at large scales. 


2 The average Mach number is M$ = b rms /C s » 0.1, where Cs is the 
sound speed, while M'^ ax a 0.3. 

3 A definition of shear and pseudo polarizations is given in section 4.2.4. 
For fluctuations with mainly perpendicular wavevectors, pseudo polarizations 
have components parallel to the mean field, while shear polarization have 
components perpendicular to the mean field. In this limit, the former are 
absent in the Reduced MHD formalism. 


along with equipartition between magnetic and kinetic energy 
and between Elsasser energies. This simulation can be classi¬ 
fied as a case of weak anisotropic turbulence in terms of the 
strength parameter^ (see table 1), although it does not have 
the properties of classical weak turbulence (Ng & Bhattachar- 
jee 1997; Galtier et al. 2000; Meyrand et al. 2014). Indeed, the 
3D spectrum has a relatively strong excitation in the parallel 
direction, resulting in a peculiar anisotropy E(k, 6) — A(6)k p , 
with an isotropic spectral index p = -2 - 3/2 in all direc¬ 
tions (corresponding to a ID spectrum with slope -3/2) and 
all the anisotropy appearing as a power anisotropy at large 
scales AW). We will not discuss the properties and the ori¬ 
gin of such a spectrum that can be found in Muller & Grap- 
pin (2005); Grappin & Muller (2010); Grappin et al. (2013)); 
what is mostly relevant for the present analysis is that run C 
has a different 3D anisotropy compared to run B, although 
in both runs energy resides mainly in perpendicular wavevec¬ 
tors. 

We will use three measures to characterize the simula¬ 
tions: II-order structure functions computed in the frame de¬ 
fined by the local mean-field (local S), II and III order struc¬ 
ture functions, respectively S and Y, computed in an abso¬ 
lute frame attached to the x-axis, which is the direction of 
the global mean-field when it is present. For all structure 
functions, computation is made calculating increments in real 
space. Local S are obtained following method I of (Cho & 
Vishniac 2000), i.e. the local field at scale t is defined as 
B f n (x) = 1/2 [b(x + t) + b(x)]. Note that the measure of 
anisotropy, i.e. the ratio S({±, 0)/S(0, t\\) of local S, is not 
unambiguously defined. In a turbulent medium fluctuations 
have a wide range of excited scales and the definition of the 
mean field depends on both scale and position. Thus, higher- 
order statistics may be introduced in the local S to a different 
extent, depending on the averaging procedure. However, the 
two-point average employed in this work is a good working 
definition, at least in simulations of homogeneous turbulence 
since it was shown to yield the same results as line average or 
volume averages, provided that the averaging scale is smaller 
than the correlation length (Matthaeus et al. 2012). The com¬ 
putation of local S is made for the whole range of available 
increments t but the average is made on a subset of grid points 
(typically N x x N y x N z = 32 3 ), which was checked to be 
a sufficient statistics for the anisotropy to converge. On the 
other hand the Ill-order structure functions Y are signed quan¬ 
tities and their computation requires large statistics to con¬ 
verge. Thus given the number of grid points N in a given 
direction, we compute Y either on a smaller range of incre¬ 
ments ({ = (1 ...N/M)dx with typically M — 4, dx is the grid 
size), or on a coarser grid of increments (t = (1 ...N/M)Mdx, 
with typically M = 4), but still on all the grid points, so the 
averages are made on a sample of > 10 7 data. The dissipation 
rate (e) entering the KHYPP equation (2) is obtained directly 
from the 3D Fourier spectra of magnetic and velocity field (b 
and u respectively), 

e = -d t E = VLkk 2 u 2 k + rfL k k 2 b 2 . (10) 

The terms appearing in Equation (2) are evaluated for four 
consecutive snapshots separated by approximately a nonlinear 
time (for the time derivative d,S we use a simple first-order 
scheme). All quantities entering the KHYPP equation are nor¬ 
malized to the dissipation rate e and then averaged over the 
four snapshots. The 3D data are finally reduced for purposes 
of representation and analysis by performing an isotropization 
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run A : SF anisotropy run B : SF anisotropy 




run C : SF anisotropy 



Fig. 1.— Anisotropy of II-order structure functions £±_{£\\) obtained by identifying the scales at which parallel and perpendicular II-order SF have the same level, 
S || = 5(^11,0) = 5(0, = S j_. Black diamonds and gray squares indicate the anisotropy with respect to the global or to the local mean-field respectively. The 

vertical bars bound inertial range scales for the global SF as determined by the Ill-order SF (see text). The two reference straight lines indicate scale-independent 
anisotropy {f\\ = dashed line) and the critical balance scale-dependent anisotropy (€\\ = , dotted line). From left to right, run A (isotropic, decaying 

turbulence), run B (anisotropic, forced, strong turbulence) and run C (anisotropic, forced, weak turbulence). Insets display the same plots compensated by ■ 


(averaging over polar and azimuthal angle 6 and 0 in spherical 
coordinates), 


v • Y\ iso 


J_ f(dYx My 

47T J \ d£ x + dly 


+ 



sin$d0d(9. 


( 11 ) 


or an axisymmetrization (averaging over the azimuthal angle 
tp in cylindrical coordinate with axis along the £ x ). 


V • FI ■ = 

T *■ I axis 



M y d£.J 


dtp. 


( 12 ) 


In the following we drop the subscripts iso and axis and even¬ 
tually mention explicitly the average procedure used for rep¬ 
resentation. 


4. RESULTS 

4.1. Anisotropy of II-order structure functions 

We measured the anisotropy of II-order SF in two frames. 
In the global frame the increments t\\ and f ± are taken parallel 
and perpendicular to a fixed direction x, which is the direction 
of the mean-field B o when it is present. In the local frame the 
parallel and perpendicular directions are relative to the scale- 
dependent mean-field direction B jj (see section 2). The mea¬ 
sure of the anisotropy is obtained by identifying the couples 
of increments (£\\, (f) at which the parallel and the perpendic¬ 
ular SF have the same value: S' u = S (C\\, 0) = S ± = 5(0, £±), 
i.e. the function €\\(f±) measures the aspect ratio of isolevels 
of the SF at different scales, also known as eddy shape. Scale- 
independent anisotropy results in a linear relation £\\ oc ( ± > 
that is, an aspect ratio £±/£\\ that does not change with scale. 
Conversely, scale-dependent anisotropy results in a deviation 
from the linear scaling £\\ oc £ p ± , with p < 1: the aspect ratio 
of SF, £±/£\\ increases at smaller scales, that is, eddies become 
more and more elongated in the parallel direction. In figure 1 
we plot (’u(Cl) for the two measures of anisotropy, local (gray 
squares) and global (black diamonds), for the three runs listed 
in table 1. The dashed line, with slope -1, is a reference for 
scale-independent anisotropy. The dotted line is a reference 
for the scale-dependent anisotropy predicted by the critical 
balance relation (Goldreich & Sridhar 1995): 


The insets display the same plots compensated by the critical 
balance anisotropy in order to better appreciate the scaling 
relation £\\(£±). 

In the local frame, all runs have a scale-dependent 
anisotropy extending to a wide range of scales. The function 
£\\{£±) has a slope flatter than 1, thus the anisotropy grows with 
decreasing scales and eddies are more and more elongated in 
the parallel direction. The scaling law actually follows the 
critical balance anisotropy Equation (13) in the range where 
S x °c £ 2 / 3 , extending for about a decade for runs A, B, and 
C in the intervals 0.008 < £± < 0.05, 0.02 < £± < 0.2, 
and 0.008 < £± < 0.08 respectively. This critical balance 
anisotropy is quite robust since in all runs the anisotropic rela¬ 
tion holds rather well even outside the above-mentioned range 
of scales where S y and S ± have a clear power-law scaling. 
Note that run C is a weak turbulence simulation, and it is 
not obvious that it should have a critical balance anisotropy 
(see Galtier et al. 2005 for an explanation based on a heuristic 
model of anisotropic turbulence). 

Consider now the anisotropy measured in the global frame 
(black diamonds in Figure 1), which will be more relevant for 
the following analysis, since it is related to the Ill-order SF 
by the KHYPP equation (2). As expected, run A is perfectly 
isotropic, the aspect ratio is unity at all scales. The anisotropy 
of strong turbulence, run B, becomes scale-independent (it 
has a slope equal to one) for scales £ ± < 0.08, approaching 
a constant ratio £\\/£x — A as 10 at small scales. For weak 
turbulence, run C, the anisotropy becomes scale-independent 
only at very small scales (£ ± < 0.01), with an aspect ratio 
A as 5 that is smaller compared to strong turbulence. The 
vertical bars in the plot bound the inertial range as identified 
from Equation (3) 4 . While run B has a scale-independent 

4 Following Equation (3), the inertial range is defined by the scales at 
which |V ■ F| is constant. We measure the slope of V ■ Y in logarithmic scales 
along the parallel and perpendicular directions and define inertial range scales 
as those ones having a slope <0.1 (see Figure 3c and Figure 4c respectively). 
This is the procedure used for run B and C. We anticipate that for run A the 
divergence of the Ill-order SF is not a nice straight horizontal line (Figure 2c), 
so in this case the inertial range is defined by the scales at which |V • F| is 
about twice larger than all other terms in the KHYPP equation (corresponding 
to a slope < 0.4). As a partial cross-check, the chosen thresholds return a 
dissipative scale that matches the scale at which Sx/(± has a local maximum, 
which is the standard procedure for the identification of the dissipative scale 
via II-oder SF. 


,2/3 


(13) 
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Fig. 2. — Run A, Bo = 0, isotropic decaying turbulence. Panel (a). Ill-order structure function, or Yaglom flux Y, projected on the {£\\, £±) plane and normalized 
by the scalar increment l. Arrows are colored according to their angle Or with respect to the radial direction (black is for 0° < Or < 5°, violet is for Or > 5°), their 
length is proportional to \Y\/€. Panel (b). Isocontour of —V • Y normalized by the dissipation rate 4e. Panel (c). Comparison of different terms appearing in the 
KHYPP equation, Equation (2), after isotropization and normalization by 4e\ the divergence of the Yaglom flux -V • Y (black solid line), the time dependent term 
d t S (dotted black line), the dissipative term 2 vd^S (red triple-dot-dashed line), and the sum of the three terms (long-dashed blue line). The thick solid horizontal 
line is a reference for 4e. The gray area in panel (a) and the white thick lines in panel (b) bound the scales at which - V ■ Y is larger than the other terms in the 
KHYPP equation, it is a rough estimate of inertial range scales. 


anisotropy that develops in the inertial range, in run C scale- 
independent anisotropy is attained only at dissipative scales. 
Thus, at inertial-range scales the anisotropy is mostly scale 
dependent in our weak turbulence simulation (run C): it fol- 
lows the critical balance scaling t\\ <x even when mea¬ 
sured in the global frame, contrary to expectations (Chen et al. 
2011 ). 

4.2. Ill-order structure functions and KHYPP equation 
4.2.1. Isotropic case 

We consider first the isotropic case (run A) for which S is 
isotropic, so we expect also to find an isotropic Ill-order SF. In 
Figure 3, panel (a), we plot the Yaglom flux Y (Ill-order SF), 
averaged along the polar angle of cylindrical coordinates with 
axis ( x = l\\, and normalized by the scalar increment I. We 
consider relatively small scales (the largest scale is l - 0.5) 
to highlight inertial range features, as will be clearer below. 
The Yaglom flux is almost radial at large scales and becomes 
remarkably radial at smaller scales £ < 0.08. The length of 
the arrows increases toward the origin, indicating that the in¬ 
tensity of the cascade increases when approaching the inertial 
range (the gray area) while keeping the same (radial) direc¬ 
tion. Note also that the arrow length is constant on circles of 
given I, meaning that Y ~ Const x l as expected for isotropic 
turbulence. Equation (6). 

In panel (b) we plot the isocontours of the divergence of 
the Yaglom flux, normalized by the dissipation rate. The iso¬ 
contours are roughly isotropic at large scales and become per¬ 
fectly isotropic at small scales (( < 0.03). In a small interval 
of scales around I as 0.01, the divergence -V • Y has a max¬ 
imum reaching the value ~ 0.9. Thus, the dissipation rate is 
approximately equal to the cascade rate, and these scales can 
be identified as the inertial range of turbulence. However, the 
divergence is not strictly a constant, as expected for the in¬ 
ertial range. Note that although the latter is very short, it is 
uniformly distributed among scales. 

Finally, in panel (c) we plot in logarithmic scales, after 
isotropization and normalization by 4e, all the terms appear¬ 
ing in the KHYPP equation, Equation (2), namely, the diver¬ 
gence of the Yaglom flux -V ■ Y (thick solid line), the dissi¬ 


pative term 2vd 2 ( S (red triple-dot-dashed line), and the time- 
dependent term d,S (dotted line). The dissipative term dom¬ 
inates at small scales, while the time-dependent term (decay) 
dominates at large scales. The cascade term, -V • Y, is larger 
than the other terms for scales 0.003 < l < 0.03 (the gray 
area in panel (a)). These two extrema can be identified with 
the injection scale and the dissipative scale, respectively. It 
is worth noting that the dissipation scale defined in this way 
coincides with the estimate based on II-order SF in Figure 1. 
From this logarithmic plot it is clear that the inertial range 
is quite small in this decaying simulation because -V ■ Y is 
much larger than other terms only in a small interval centered 
at ( as 0.01, where it is roughly horizontal and equal to 0.9 (it 
should be equal to one in an ideally infinite inertial range). 

In the same plot we also traced, as a blue long-dashed line, 
the sum of the three terms just discussed, that should amount 
at all scales to the dissipation rate 4e (thick solid horizon¬ 
tal line) for good energy conservation. The sum is an almost 
horizontal straight line, only a factor 1.2 higher than the dissi¬ 
pation rate for scales t > 0.005. This confirms that the statis¬ 
tics is large enough to ensure convergence and that the con¬ 
servation of energy holds with sufficient accuracy except at 
very small scales where a small numerical dissipation proba¬ 
bly kicks in: the injection at large scale (including the decay), 
the cascade in the inertial range, and the dissipation at small 
scales all occur at the same rate. 

To summarize, the analysis of the isotropic turbulence con¬ 
firms the theoretical expectations: inside the inertial range 
the Yaglom flux is radially directed and its magnitude scales 
linearly with the scalar increment £. The divergence of the 
Yaglom flux is approximately constant in the inertial range, 
and it is uniformly distributed among scales (isotropic). Note, 
however, that the extent of the inertial range is very limited 
due to the relatively small Reynolds number that prevents the 
formation of a large range of scales where the divergence 
of the Yaglom flux is the dominant term. With the current 
resolution (1024 3 ) it is at best half a decade, indicating that 
this is the minimal resolution required for studies of decaying 
turbulence (although hyperviscosity would probably alleviate 
the problem). This is relevant for solar wind studies (Dong 
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Fig. 3.— Run B. Bo = 5, strong turbulence. Panel (a). Ill-order structure function, or Yaglom flux Y, projected on the (t\. ( ± ) plane and normalized by the 
perpendicular increment (±. Arrows are colored according to their angle 0± with respect to the perpendicular direction direction (black is for 0 ° < 6± < 5°), their 
length is proportional to |Y|/f_L- Panel (b). Isocontour of -V ■ Y normalized by the dissipation rate 4e. Panel (c). Cuts of -V ■ Y in directions parallel (dashed 
line) and perpendicular (solid line) to the mean-field Bo. The thick solid horizontal line is a reference for 4e. The gray area in panel (a) and the white thick lines 
in panel (b) bound the scales at which -V ■ Y is larger than the other terms in the KHYPP equation, it is a rough estimate of inertial range scales. 


et al. 2014), in which expansion induces an additional decay 
of magnetic and kinetic energy on top of the decay due to the 
turbulent dissipation. 

4.2.2. Anisotropic case, strong turbulence 

Let’s now turn to the anisotropic run B, which is a simula¬ 
tion of strong turbulence with guide field and forced at large 
scales on components perpendicular Bq. In Figure 3, panel (a) 
we plot again the Ill-order SF, i.e. the Yaglom flux Y, aver¬ 
aged over the polar angle in cylindrical coordinates with axis 
along increments parallel to the mean magnetic field. At vari¬ 
ance with the isotropic case, the arrows are colored according 
to the angle 9± formed with respect to the perpendicular di¬ 
rection (black color stands for 0 ± < 5°), and their length is 
normalized to the perpendicular increment £ ± = {t 1 + f 2 ) 1 ' 2 . 
The Yaglom flux is remarkably vertical in the inertial range 
(the gray area), and it is proportional to the perpendicular in¬ 
crements Y oc f ± (the arrow length is uniform in the whole 
inertial range, after normalization). This suggest that turbu¬ 
lence is undergoing a purely 2D cascade. Equation (8). 

In panel (b) the normalized divergence of the Yaglom flux is 
constant over a large interval of parallel and perpendicular in¬ 
crements, with a value close to 1 (the light-green area at level 
0.9). However it is not uniformly distributed among scales, 
the isocontour of constant divergence extends to smaller per¬ 
pendicular scales, suggesting that the cascade is not removing 
positive correlation from the parallel direction. This is consis¬ 
tent with Figure lb that shows a clear anisotropy of S in favor 
of a two-dimensionalization in the perpendicular plane. 

This can be better appreciated in panel (c), where we plot 
cuts along the parallel and perpendicular directions of the di¬ 
vergence of the Yaglom flux in logarithmic scales. There is a 
very nice constant divergence in the perpendicular direction, 
with a value close to the dissipation rate (the horizontal thick 
solid line at level 1), covering about one decade in the range 
0.01 < {± < 0.1. In the parallel direction, an approximate in¬ 
ertial range is also found, having a smaller extent and shifted 
to larger scales 0.1 < t\\ < 0.6. 

We recall that from -V • Y (panels (b) and (c)) one can 
identify the location of the inertial range and its distribution 
among parallel and perpendicular scales. On the other hand, 
the simple dependence of the Yaglom flux allows us to iden¬ 


tify the cascade as a 2D cascade with Y = -2et ± We antic¬ 
ipate that Y + = ^Az ± |Az + | 2 ^ have only perpendicular com¬ 
ponents because there is a dominance of shear-Alfven polar¬ 
izations. Indeed, these polarizations have Az ± lying in the 
perpendicular plane and, as a consequence, the cascade is 2D. 
This simple picture of 2D cascade does not hold anymore as 
soon as pseudo polarizations, which have out-of-plane com¬ 
ponents, are energetically important. 

4.2.3. Anisotropic case, weak turbulence 

We finally consider the case of weak turbulence with guide 
field (run C), which is forced isotropically at small wavevec- 
tors 1 < |&| < 2 by imposing at all times (freezing) the cor¬ 
responding Fourier modes of the fluctuating fields B , u (note 
that their x,y,z components have also equal energy). In fig.4, 
panel (a) one can see immediately that the Yaglom flux is 
oblique, with an angle 8 ± that changes with scales (we nor¬ 
malize arrow length by the perpendicular increment |F|/£_l as 
in Figure 3). This means that the Ill-order SF has a parallel 
component and a non-negligible dependence on the parallel 
increments, thus contributing to the cascade rate through the 
divergence of the Yaglom flux. Such contribution seems to 
decrease at small parallel scales, where the Yaglom flux be¬ 
comes vertical, hinting to a milder two-dimensionalization of 
this weak turbulence cascade. 

Isocontours of -V • Y are plotted in panel (b), with the 
usual normalization by 4e. The inertial range can be iden¬ 
tified with the light-orange area at level 1.1, extending to a 
wide range of parallel and perpendicular scales. Its distribu¬ 
tion is non-uniform and more complex than that of strong tur¬ 
bulence (Figure 3b), reflecting the dependence of the Ill-order 
SF from t ± and €\\. 

Panel (c) shows cuts of the divergence of the Yaglom flux 
along the parallel and perpendicular directions. Although the 
divergence is not exactly constant, the inertial range covers 
more than one decade in the perpendicular direction, 0.005 < 
( ± < 0.1, yielding a cascade rate that is slightly higher than 
the dissipation rate 4e. On the other hand, the cut in the paral¬ 
lel direction is less flat, making more questionable the identi¬ 
fication of an inertial range in this direction. The approximate 
parallel inertial range is shorter and located at larger scales, 
0.03 < 4 < 0.3. 
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Fig. 4.— Run C, Bo = 5, weak turbulence. Same as in Figure 3. In panel (a) the arrow colors, from black to purple, indicate € [0“, 90"| binned in intervals 
of 5°. 


4.2.4. r/te KHYPP equation for pseudo and shear Alfven 
polarizations 

Decomposing fluctuations in pseudo and shear Alfven po¬ 
larizations proves useful to analyze in some more detail the 
relation between the Yaglom flux and the cascade rate in run 
C, which is a simulation of incompressible MHD. The decom¬ 
position is made in Fourier space, where pseudo Alfven po¬ 
larizations and shear Alfven polarizations are oriented along 
the unitary vectors (e.g. Maron & Goldreich 2001), 

kxB 0 B 0 -(k■ B 0 )k 

* sh [1 - (k ■ Bo) 2 ] 1 ' 2 ’ * ps [1 -(k-Bp) 2 ] 1 ' 2 ' ( 

This decomposition is completely equivalent to the decompo¬ 
sition into toroidal and poloidal components of the magnetic 
fluctuations. In incompressible MHD the same decomposi¬ 
tion applies to the velocity field, which is also solenoidal. 
Shear Alfven polarizations are the proper Alfven modes in 
full MHD. Their component is perpendicular to both the mean 
field B o and the wavevector k, being incompressible. The 
pseudo Alfven polarizations are the incompressible limit of 
slow modes in MHD. Their component lies in the plane iden¬ 
tified by the mean field and the wavevector, and it is again 
perpendicular to the wavevector. For fluctuations with strong 
anisotropic spectra (k± » k||), shear polarizations have 
wavevectors and components lying in the plane perpendicu¬ 
lar to Bp, thus they represent the 2D modes in the slab-plus- 
2D decomposition introduced by Matthaeus et al. (1990). In¬ 
stead, pseudo polarizations have wavevectors in the plane per¬ 
pendicular to Bo but components along Bo- This polarization, 
which is absent in Reduced MHD, is instead present in 2.5D 
configurations with out-of-plane mean field, and should not 
be confused with the slab component that has wavevectors 
parallel to Bq and fluctuations perpendicular to it (and is thus 
included in Reduced MHD). 

After decomposing fluctuations in pseudo and shear Alfven 
polarizations, we go back to the real space and compute sep¬ 
arately all the contributions to the KHYPP equation. Equa¬ 
tion (2). Note that for strictly parallel wavevectors the pseudo¬ 
shear decomposition degenerates, so we remove such modes 
(the slab component) in the following analysis to avoid arbi¬ 
trary partition of energy into the shear and pseudo polariza¬ 
tions. Using Parseval’s theorem (Az* ■ Az* ) = 0, the decom¬ 


posed KHYPP equation can be written as: 


(S sh + S ps) — 2 vV^(S s h + S ps ) + 4e — 

(15) 

~ * {^AZsh\AZsh\ "t" AZsh\AZps\ ^ 

(16) 

— ‘ (^AZps\AZsh\ + ^Zps\AZps\ ^ 

(17) 

— ^£ * ^2A Zsh(AZsh * AZps) 2Az ps (AZsh 

■A z P s)), (18) 


where we summed the ± species and dropped ± super¬ 
scripts, i.e. S sh = l/2<|Az+ ft | 2 + |Az~ | 2 ), Az sft |Azp S | 2 = 

l/2[Az+JAz“J 2 + Az“JAZpJ 2 ], etc ... 

The Ill-order SF on the rhs is split into three lines contain¬ 
ing respectively (1) the strain of the shear polarizations on the 
shear and pseudo energies, (2) the strain of the pseudo polar¬ 
izations on the shear and pseudo energies, and (3) the mixed 
terms accounting for the exchange of energy during the cas¬ 
cade between the shear and pseudo polarizations. We antici¬ 
pate that the mixed terms are negligible in the inertial range, 
thus we split the above equation into a system of equations for 
S ps and S s h, with their own cascade rate e v /, and e ps , 

@lS sh ~ 2vV f S s h + 4 6 s h — ~^t ’ (AZsaIAZsaI + Az ps \AZsh \~^ 5 

(19) 

2vV^5 ps + 4 € ps — ~• ^AZs/i|AZps| - 3" AZps|Az p s|-^ . 

( 20 ) 

Note that each equation contains only quadratic terms of 
a given Alfven polarization (shear or pseudo) even in the 
Ill-order SF. On the rhs of Equations (19)-(20), the sec¬ 
ond terms represent an active contribution to the cascade of 
pseudo Alfven polarizations: if it is vanishing or negligible 
the pseudo Alfven polarizations can be said passive. 

This separation will appear justified by the analysis of each 
term, which is presented in Figure 5 in the same format of 
figs. 3, 4, i.e. with the the Yaglom flux normalized by the per¬ 
pendicular increment, Y/£ ± , and the divergence normalized 
by the cascade rate, -V- Y /4e. We consider separately the two 
contributions in the rhs terms of the equation for the pseudo 
energies. Equation (20). In the first column we have the Ill- 
order SF accounting for the strain of shear polarizations acting 
on the pseudo energy, Y sp = ^Az s ft|Az ps | 2 ^ (the first term in 
the rhs). The Yaglom flux (upper panel) is perpendicular, as 
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Fig. 5.— Run C. Decomposition of the KHYPP equation for shear and pseudo energies. As in figs. 3, 4, Y is normalized by £j_, -V • Y is normalized by 4e. 
In the first three columns we analyze the KHYPP equation for pseudo energy. Equation (20), by plotting the contribution to the cascade appearing on the rhs of 
Equation (20). The first term, i.e. the strain of shear polarizations on the pseudo energy (column 1); the second term, i.e. the strain of pseudo polarizations on 
pseudo energy (column 2); and their sum (column 3). In column 4 we consider the cascade for shear energy. Equation (19), without separating the pseudo and 
shear contribution. In the bottom panel, fourth column, we also plot with grey color the contribution of the mixed terms. Equation (18). 


in run B, since the vectorial part of Y sp is made of only shear 
polarizations, which have only perpendicular components. Its 
magnitude depends mostly on £ ± as can be guessed by the 
length of the arrows, with a small dependence on £\\ at large 
parallel scales. This dependence is better appreciated on the 
map of -V • Y sp (middle panel) that has slightly inclined iso¬ 
contours instead of horizontal isocontours. The inertial range 
identified with the green area at level * 0.4 is non uniformly 
distributed in the (£\\, £±) plane, occupying preferentially per¬ 
pendicular scales £ ± < 0.05 independently of parallel scale. 
In the bottom panel one can better locate the inertial range 
which is much more extended toward smaller scales in the 
perpendicular than in the parallel direction. 

In the second column, we analyze the Ill-order SF account¬ 
ing for the strain of pseudo polarizations on the pseudo en¬ 
ergy, y pp = ^Azp S |Azps| 2 ), the second term in the rhs of Equa¬ 
tion (20). Note that the Yaglom flux is now horizontal (ex¬ 
cepts at small parallel scales £\\ < 0.03), since for a spectrum 
(or II-order SF) with energy residing mainly in perpendicular 
wavevectors, pseudo polarizations have a dominant parallel 
component. The contribution to the divergence (middle panel) 
is almost complementary to the previous term, being larger at 
large parallel and perpendicular scales while it is negligible 


in the inertial range. This is better seen in the bottom panel, 
in which it is also shown clearly that this Ill-order SF does 
not form an inertial range on its own, but it is responsible for 
a non-constant energy flux from large to small scales in both 
parallel and perpendicular directions. One would be tempted 
to model the two terms, y sp and y pp , as 2D+1D components 
respectively, according to the orientation of the Yaglom flux. 
However, such a decomposition does not work since the can¬ 
didate for the ID component (Y pp ) has a strong dependence 
on £ ± (second column, top panel), and its divergence does not 
yield any identifiable inertial range (second column, middle 
and bottom panels). 

Finally in the third column we analyze the sum of 
the two terms just discussed, Y* p = Y pp + Y sp = 
((Azsi, + Azp S )|Az i , s | 2 ), the whole rhs of Equation (20). The 
Yaglom flux is now oblique, with dependence on both £ ± 
and £||, the latter being mostly inherited from the strain of 
pseudo polarizations (Y pp ), indicating a smaller degree of 
two-dimensionalization for turbulence of pseudo Alfven po¬ 
larizations, compared to the strong turbulence of run B. In the 
middle panel the inertial range extends to parallel and per¬ 
pendicular scales in a way somewhat similar to run B, al- 
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though the Yaglom flux is quite different. Comparing the 
bottom panels in columns one, two, and three, one can see 
that the strain of both shear and pseudo polarizations con¬ 
tribute to the formation of the wide inertial range in the per¬ 
pendicular direction 0.005 < £ ± < 0.1, as well as to the 
shorter inertial range found at larger scales in the parallel di¬ 
rection 0.03 < £± < 0.2 (their ranges coincide with those one 
obtained without separating shear and pseudo polarizations). 
More precisely, shear polarizations contribute with a constant 
cascade at small scales (£ < 0.04), while pseudo polarizations 
control the injection of energy into the cascade at large scales 

« 0.1). 

We do not repeat the analysis of the strain of shear and 
pseudo polarizations on shear energy. Equation (19), since it 
yields qualitatively similar results. The only noticeable differ¬ 
ence is that the Yaglom flux Y ps = ^Az p *|Az sft | 2 ^ is radially 
directed instead of being horizontal. We plot in the fourth 
column the whole contribution to the cascade of shear energy, 
y* s = ^( Az s h + Az P s)|Az s / ! | 2 ^. Comparing the third and fourth 
columns, one can see that (i) shear and pseudo energy cas¬ 
cade qualitatively in a similar way (top panels), although the 
Yaglom flux for shear polarizations, Y* s , is more perpendic¬ 
ular (more 2D); (ii) the cascade rate is slightly stronger for 
shear energy than for pseudo energy (middle panels); and (iii) 
a clear inertial range is found for both shear and pseudo en¬ 
ergies (bottom panel) following the decomposition of Equa¬ 
tions (19)-(20). 

In the fourth column, bottom panel, we also plot the sum of 
the mixing terms. Equation (18), along the parallel and per¬ 
pendicular increments (gray lines). They are negligible com¬ 
pared to all other terms in the inertial range, indicating that 
the shear and pseudo polarizations cascade without exchang¬ 
ing their energy, as already pointed out by Maron & Goldre- 
ich (2001) for decaying simulations. However, at large scales 
£ ~ 0.2, they are of the same order of the term accounting for 
the strain exerted by pseudo polarizations (II column, bottom 
panel in Figure 5), suggesting an exchange of energy between 
the shear and pseudo polarizations. 

To summarize, in run C the freezing of large scales (forc¬ 
ing) causes an exchange of energy between pseudo and shear 
polarization at large scales. The two polarizations cascade 
in a similar manner (same rate, same inertial range) under 
the action (strain) of both the shear and pseudo polarizations. 
The flux toward small scales that is triggered by pseudo po¬ 
larizations is smaller in magnitude and not constant (it is not 
a proper inertial range on its own), but it affects considerably 
the total cascade rate at large inertial-range scales (£ ~ 0.1). 
On the other hand, the shear polarizations control the constant 
energy flux at small inertial-range scales (£ < 0.04). Thus ac¬ 
cording to the separation made in Equations (19)-(20), the two 
cascades are not totally independent. Indeed, whatever the po¬ 
larization considered (S ps or S v /,), the cascade is triggered by 
pseudo polarizations at the interface between injection scales 
and inertial range scales, while it is sustained by shear po¬ 
larizations at smaller inertial range scales. The presence and 
importance of the pseudo polarizations is revealed by the Ya¬ 
glom flux, the more pseudo polarizations are energetic, the 
more Y becomes oblique. 

5. ANISOTROPY OF THE III-ORDER SF AND SOLAR 
WIND MEASUREMENTS 

In measuring the cascade rate of solar wind turbulence 
through Ill-order SF, one needs to make assumptions on the 


unknown anisotropy of Y. Despite the limitations due to the 
imposed symmetries, this brings the advantage of increasing 
the statistic, on the one hand, and of avoiding the computa¬ 
tion of the divergence of Y on the other hand. Two common 
assumptions usually employed to reduce in-situ data will be 
tested against our anisotropic DNS (run B and C) in order 
to highlight possible systematic errors on the measure of the 
cascade rate in the solar wind. 

The simplest assumption is that of isotropic turbulence. 
Equation (6), whose expression is repeated here for conve¬ 
nience. 


--III 

4 £ ’ 


( 21 ) 


and involves projecting Y along the direction of increments 
(Ye = Y ■ (If). As discussed in section 2, the cascade rate 
is directly obtained from the Ill-order SF without any need 
to compute derivatives. This “isotropic” method has been 
mainly employed for fast polar wind (Sorriso-Valvo et al. 
2007; Marino et al. 2008, 2012) on Ulysses measurements. 

Admitting some form of anisotropy, the minimal assump¬ 
tion is that of axisymmetry. To avoid computing derivatives 
along the two increments, one again resorts on simplified 
geometrical models, as the 2D+1D turbulence employed by 
Macbride et al. (2008); Stawarz et al. (2009, 2010) on WIND 
and ACE data in the ecliptic solar wind. This method assumes 
that the Ill-order SF has parallel and perpendicular compo¬ 
nents that depend only on the parallel and perpendicular incre¬ 
ments, respectively (the cascade is independent in the two di¬ 
rections). The total cascade rate is obtained by combining the 
two independent equations for (isotropic) 2D-perpendicular 
and for ID-parallel cascades (Equation (8) and Equation (9) 
respectively). Such a “hybrid” cascade rate reads 


^ = e .» t6 » = e . + ^ = _(I3 + >A). (22) 

The isotropic and hybrid cascades are obtained by first com¬ 
puting the 3D Ill-order SF, and then by applying the corre¬ 
sponding average (isotropy or axisymmetry). Similarly, the 
“true” cascade rate is obtained from Equation (12): 


^true 


n 

4 2ttJ \d£ x d£ y 



(23) 


We now compare in Figure 6 the true cascade rate (thick 
solid line) with the isotropic (thin solid line) and the hybrid 
(thin dashed line) cascade rate for run B and C (top and bot¬ 
tom panels, respectively). All cascade rates are normalized 
by the dissipation rate e obtained directly from spectra. Equa¬ 
tion (10). Each panel represents a cut in the (£\\, £ ± ) plane 
taken along a fixed direction that forms an angle 9 with the 
mean-field direction. 

The isotropic cascade rate is strongly angle dependent, in¬ 
dependently of the run considered. It returns increasing cas¬ 
cade rates for increasing angles, as can be expected since both 
runs B and C have a dominant perpendicular cascade. It over¬ 
estimates the true cascade rate at large angles (8 > 70°) while 
it underestimates the true cascade rate by a factor « 2 even at 
8 - 45°. The hybrid method instead performs very well: it 
does not vary with angle 8 and yields the correct cascade rate 
at oblique angles 8 > 20". 

In the nearly parallel direction, both the hybrid and 
isotropic methods strongly underestimate the true cascade rate 
(by a factor > 10) or they completely fail in yielding any 



11 


run B : e(true)/«, e(iso)/e, e(hybr)/« 



Fig. 6. — Comparison of the true cascade rate, Equation (23) (solid thick line) with the isotropic cascade rate. Equation (21) (thin solid line), and the hybrid 
cascade rate. Equation (22) (dashed black line). The red dashed line is the ID (parallel) cascade rate obtained within the hybrid method. From left to right, panels 
refer to increasing angles 8 between the sampling direction l and the mean-field Bo. 


cascade rate. Indeed, neither of the two methods is able to 
account for the dependence of the parallel component of the 
Ill-order SF on the perpendicular increments, Y\\{£±), which is 
instead fundamental in our strongly anisotropic runs (e.g. Fig¬ 
ure 2a). Indeed, the ID cascade entering the hybrid method 
is basically negligible at all angles (red dashed line in fig 6). 
This implies that the cascade is with a good approximation a 
2D cascade in both our anisotropic runs. On the other hand, 
a linear scaling for the ID cascade is found in the solar wind 
(e.g. Macbride et al. 2008), suggesting that the Ill-order SF of 
solar wind turbulence is less anisotropic than that one of runs 
B and C. 

6. SUMMARY AND DISCUSSION 

We studied the anisotropy of MHD turbulence with or with¬ 
out guide field (Z?o) by means of structure functions (SF) of II 
and III order, in three direct numerical simulations of MHD 
turbulence at moderate and high resolution (see table 1). We 
consider isotropic decaying turbulence (run A), forced strong 
turbulence with mean-field (run B), and forced weak turbu¬ 
lence with mean-field (run C). 

We used II-order SF ( S ) to characterize the anisotropy with 
respect to the local mean-field and to the global mean-field. 
The anisotropy with respect to the local mean-field is scale 
dependent in all runs (with or without mean-field) and follows 
the critical balance prediction (\\ oc l[ . This confirms previ¬ 
ous theoretical and numerical findings (Goldreich & Sridhar 
1995; Cho & Vishniac 2000; Maron & Goldreich 2001), in 
the local frame the anisotropy grows at smaller and smaller 
scales. When instead we used as a reference the global mean- 
field, we found isotropy for MHD turbulence without mean- 
field (run A) and strong two-dimensionalization for turbu¬ 
lence with mean field, with run B being more anisotropic than 
run C (the small-scale aspect ratio is 10 and » 5 respec¬ 
tively). Surprisingly, for run C, we found a scale-dependent 
anisotropy consistent with critical balance even when S is 
computed in the frame attached to the mean field Bo, at vari¬ 
ance with previous analysis in DNS (Chen et al. 2011) and 
in solar wind measurements (Tessein et al. 2009) which re¬ 
ported scale-independent anisotropy. A possible explanation 
is that run C has a strong magnetic field (b rms /Bo = 1/5) and 
so the local and global frames almost coincide. In favor of 
this explanation, the anisotropy becomes scale independent 
at smaller scales. However, there are other “non-standard” 
aspects of run C that could be related to scale-dependent 


anisotropy in the global frame: the special 3D spectrum that 
has an isotropic spectral slope but anisotropic energy levels 
in parallel and perpendicular direction (Grappin & Muller 
2010), the strong excitation of pseudo Alfven polarizations 
(see below), and the isotropic forcing (freezing) that main¬ 
tains a magnetic excess at large scales (Grappin et al. 2013). 
These three aspects could be all related to one another and are 
under investigation. 

We then analyzed the vectorial Ill-order SF (or Yaglom 
flux, Y), which is related to S computed in the global frame by 
a generalization to incompressible MHD (Politano & Pouquet 
1998) of the Von Karman-Howart-Yaglom relation in hydro¬ 
dynamics (KHYPP equation in the following). The Ill-order 
SF, Y, is expected to be anisotropic for MHD turbulence with 
mean-field, but it has never been reported in DNS or exper¬ 
iments (except for rotating hydrodynamic turbulence experi¬ 
ments, Lamriben et al. 2011). Despite some anisotropic mod¬ 
els exist (Podesta et al. 2007; Galtier 2009, 2012), they have 
never been verified in DNS. 

We tested the KHYPP equation in isotropic decaying tur¬ 
bulence (run A) and found that the Yaglom flux is almost 
purely radial, its divergence is isotropic, and the maximum 
of -V • Y is of the order of the dissipation rate (only a factor 
0.8 smaller), confirming theoretical expectations (Politano & 
Pouquet 1998). However, the condition defining the inertial 
range, -V ■ Y = const, is only fairly satisfied in our decay¬ 
ing turbulence that has an inertial range of only half a decade 
despite the high resolution (1024 3 ). More interestingly, we 
showed that the familiar concept of Kolmogorov cascade in 
Fourier space also applies in the real space with SF. Note that 
no assumption on locality of nonlinear interactions is needed 
to obtain the KHYPP equation. Despite this, we found that 
injection, cascade, and dissipation occur at the same rate and 
that the cascade is achieved via a fairly constant nonlinear 
transfer in the inertial range. In particular, by comparing the 
divergence of the Yaglom flux with the other terms appearing 
in the KHYPP equation (dissipation term and decaying term), 
one can determine the dissipation scale and the injection scale 
of turbulence. This has important applications in solar wind 
turbulence, since, due to expansion, it is an intrinsically de¬ 
caying turbulence for which the injection scale is not easily 
identified (Hellinger et al. 2013; Dong et al. 2014). 

For strong turbulence with mean-field, we found that Y 
has only perpendicular components and depends only on per¬ 
pendicular increments; this is a 2D turbulence as defined in 
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Equation (7). The corresponding inertial range extends to 
small perpendicular scales, while it is shorter and located at 
larger scales in the parallel direction: the anisotropy of Y 
emerges as a non-uniform distribution of the scales having 
-V • Y = const. A similar non-uniform distribution is found 
in the forced, weak turbulence case (run C), with an impor¬ 
tant difference concerning the orientation of Y. The latter is 
no more strictly perpendicular, but oblique. In particular, the 
Yaglom flux is almost radial at large parallel scales in the iner¬ 
tial range and becomes more and more perpendicular at small 
parallel scales. We interpret this behavior as a weaker two- 
dimensionalization of turbulence, since the oblique orienta¬ 
tion of Y implies a dependence on parallel and perpendicular 
scales, in contrast to the purely 2D case in which Y = Yj_(( ± ). 

By decomposing fluctuations into shear Alfven and pseudo 
Alfven polarizations in run C we also found that energetically 
important pseudo polarizations are responsible for the non¬ 
perpendicular Ill-order SF. The pseudo/shear decomposition 
allows us to prove that the two polarizations cascade indepen¬ 
dently, that is, without exchanging their energy, in the iner¬ 
tial range. One can thus write separate KHYPP equations for 
pseudo and shear energies. However, strictly speaking, the 
two cascades are not independent, since in each of them the 
strain of pseudo and shear polarizations enters the expression 
for the Yaglom flux. In particular the strain exerted by pseudo 
polarizations is associated with a non-constant cascade of en¬ 
ergy at the large inertial-range scales; thus pseudo polariza¬ 
tions control the injection of energy in run C (although be¬ 
coming passive at smaller inertial-range scales, e.g. Maron 
& Goldreich 2001; Cho & Lazarian 2003). In Grappin et al. 
(2013), it is shown how the equipartition of kinetic and mag¬ 
netic energy at the large frozen scales is associated to the 
weaker anisotropy of the 3D spectrum of run C, having the 
property of a 3D Iroshnikov-Kraichnan turbulence. Here we 
argue that kinetic/magnetic equipartition at the largest scales 
induces an exchange of energy between pseudo and shear po¬ 
larizations, thus allowing pseudo polarizations to control the 
cascade at the interface between injection and inertial-range 
scales. 

Finally, on the anisotropic runs B and C we tested two meth¬ 
ods that are commonly applied to obtain the cascade rate in 
the solar wind. The methods both rely on assumptions of the 
anisotropy of the Ill-order SF, allowing one to obtain the cas¬ 
cade rate directly from Y (i.e. without computing the diver¬ 
gence). The first method assumes that Y is isotropic. The sec¬ 
ond method (hybrid method) assumes axisymmetry and mod¬ 
els the anisotropy as a geometric superposition of a ID (paral¬ 
lel) cascade and an isotropic 2D (perpendicular) cascade. We 
found that the isotropic method is strongly angle dependent, 
yielding correct cascade rates only at large angles between 
the direction of increments and the magnetic field, translat¬ 
ing to angles between the magnetic field direction and the so¬ 
lar wind direction Qbv <; 70" (needless to say, the isotropic 
method works very well for run A). For fast polar wind, em¬ 
anating from stable coronal holes, turbulence is expected to 
have a strong anisotropy (Verdini et al. 2012; Perez & Chan- 
dran 2013), comparable to our run B and C. We thus suggest 
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that the angle dependence of the isotropic method is at the ori¬ 
gin of the relatively small number of intervals (13%) in which 
the linear scaling Y{ oc { is found in polar solar wind in-situ 
data (Marino et al. 2012). In contrast, the hybrid method per¬ 
forms very well on our anisotropic runs, yielding an angle- 
independent cascade rate for 8bv 20°. Both methods fails 
to obtain any cascade rate when increments are parallel to the 
mean-field. This is due to the dependence F||(Fl) that is ex¬ 
cluded in both assumptions but present in all our runs. This, 
together with the good performance of the hybrid method, in¬ 
dicates that the anisotropic runs B and C are fairly well de¬ 
scribed by a 2D cascade model, despite their different degree 
of anisotropy. Note that in the solar wind a ID cascade, ac¬ 
tually a linear scaling Y\\ oc t\\, is indeed measured with the 
hybrid method, indicating that solar wind turbulence in the 
ecliptic has a different and probably much weaker anisotropy 
compared to our DNS. 

We conclude by noting that the ratio b rms /Bo « 1/5 em¬ 
ployed in our simulation is lower than the actual value found 
in the solar wind, which is « 1/2 at scale of few hours. The 
small value of fluctuations’ amplitude was chosen to highlight 
the effects of anisotropy and makes the comparison of our re¬ 
sults with solar wind turbulence at large scales a bit weak, 
being more meaningful for scales shorter than a minute. In a 
future work we plan to apply a similar analysis to simulations 
with larger ratio b rms /Bo also dropping the assumption of ax¬ 
isymmetry. In fact. Cluster observations (Narita et al. 2010) 
showed that the solar wind turbulence is not axis-symmetric 
(but see Turner et al. 2011 for an explanation in term of an 
observational bias). DNS of MHD turbulence in the frame¬ 
work of the Expanding Box Model (Grappin et al. 1993; Grap¬ 
pin & Velli 1996) also show that expansion causes non axis- 
symmetric 3D spectra (Dong et al. 2014). Note that expansion 
is also responsible for a differential decay of fluctuations, thus 
cross helicity and magnetic/kinetic imbalance of fluctuations 
also enter the KHYPP equation in the expanding solar wind 
(Hellinger et al. 2013; Gogoberidze et al. 2013). Another in¬ 
teresting direction is to explore the role of velocity shears on 
the KHYPP equation (Wan et al. 2009, 2010), which have 
been invoked to be the main driver for the turbulence cascade 
in both the ecliptic (Roberts et al. 1992) and the polar solar 
wind (Marino et al. 2012). Again, numerical simulations with 
the Expanding Box Model seem to be a very promising tool to 
understand the combined effect of shear and expansion (e.g. 
Roberts & Ghosh 1999) in shaping the cascade of solar wind 
turbulence. 
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